Comparison

Introduction

The results below shows the comparison of several normalization and imputation methods. We decided to use 4 imputation methods:

1. Shifted Distribution based on the data;
2. Imputation using the Random Forest algorithm;
3. Imputation using the k-Nearest Neighbors algorithm;
4. Imputation using the ludovic method.

and 8 normalization methods:

1. Z-Score Normalization
2. Min-Max Normalization
3. Median Scaling
4. MAD (Median Absolute Deviation) Scaling
5. Linear Regression Normalization
6. VSN (Variance Stabilization Normalization)
7. EigenMS

We selected the following features to train the machine learning models:

• Sequence coverage - percentage of the sequence that is covered by the identified peptides of the first protein sequence contained in the group;
• Sequence length - the length of the leading protein sequence contained in the group;
• Mol weigth kDa - molecular weight of the leading protein sequence contained in the protein group;
• Qvalue - the ratio of reverse to forward protein groups;
• Score - Protein score which is derived from peptide posterior error probabilities
• MS/MS Count - the number of tandem MS scans that were made based on this full scan (e.g. a top 10 method selects the top 10 most intense ions in the scan and fragments those).
• LFQ values for TOTALS (KO 22/23/24 and WT 22/23/24) - target
• Based on the column names of tje LFQ values, we created 2 more features:
  1. celltype - Knockout (KO) or Wild Type (WT) - for kNN encoded as (1:KO, 2:WT)
  2. rep - replicates (22, 23 or 24) - for kNN encoded as (1:22, 2:23, 3:24)

We used the following parameters for the K-Nearest Neighbors algorithm (bold row):

k	RMSE	Rsquared	MAE
5	560141411	0.5764903	83499878
7	586141572	0.5110148	89345187
9	615649021	0.4549313	93976048

For the Random Forest algorithm we used the Grid Search method to find the best hyperparameters. The table belows shows the output of the top 10 models.

GRID SEARCH FOR RANDOM FOREST MODEL

	num.trees	mtry	min.node.size	replace	sample.fraction	rmse	perc_gain
1	100	4	5	false	0.8	36184321.6466	2.8631
2	500	4	5	false	0.8	367032968.8757	1.4713
3	300	4	5	false	0.8	367963900.30284	1.2214
4	500	4	1	true	0.8	370084110.1524	0.6523
5	500	4	10	false	0.8	370380581.2619	0.5727
6	300	4	1	true	0.63	371446664.9562	0.2865
7	100	4	5	false	0.8	372131148.9947	0.1027
8	300	4	10	false	0.8	372169999.3198	0.0923
9	300	4	4	true	0.63	372767670.8139	-0.0681
10	300	4	3	true	0.8	373586116.5828	-0.2878

Base on the output of the grid search, we can see that there is better model than default and it is random forest model with the following hyperparameters:

• num.trees = 100
• mtry = 4
• min.node.size = 5
• replace = FALSE
• sample.fraction = 0.80

We train Linear Regression Model which, after the stepwise selection, has the following formula:

\[LFQvalue = 8631815 \cdot no.proteins -18172773 \cdot peptides + 1794237 \cdot seq.coverage +\\ + 807933 \cdot mol.weight.kDa + 430143 \cdot score + 1490175 \cdot ms.count\]

Where:

• no.proteins - Number of proteins contained within the group.
• peptides - The total number of peptide sequences associated with the protein group (i.e. for all the proteins in the group).
• seq.coverage - Percentage of the sequence that is covered by the identified peptides of the first protein sequence contained in the group.
• mol.wieght.kDa - Molecular weight of the leading protein sequence contained in the protein group.
• score - Andromeda score for the best associated MS/MS spectrum.
• ms.count - The number of MS/MS spectra recorded in this raw file.

We train and compare 4 machine learning models and based on the performance (table below), we decided to use 2 of them:

1. Random Forest
2. kNN

Model	NRMSE	MSE	RMSE	MAE	R2
RF	0.5724	0.8044	0.8969	0.6663	0.7615
kNN	0.6070	1.0182	1.0091	0.7590	0.6701
LM	0.8737	3.8825	1.9704	1.6091	0.2623
PCA	1	6.8863	2.6242	2.2380	0

For the Ludovic method first we had to assign the missing values using the protti::assign_missingness() fxn, which added the comparison column containing the comparison name for the specific treatment/reference pair and the missingness column with 4 types of missingness:

• "complete": No missing values for every replicate of this reference/treatment pair for the specific grouping variable.
• "MNAR": Missing not at random. All replicates of either the reference or treatment condition have missing values for the specific grouping variable.
• "MAR": Missing at random. At least n-1 replicates have missing values for the reference/treatment pair for the specific grouping variable.
• NA: The comparison is not complete enough to fall into any other category. It will not be imputed if imputation is performed.

Then with additional column we could enter our data into protti::impute() fxn, which performs the Ludovic imputation method based on the missingness column using the following instructions:

• "MNAR" missingness is sampled from a normal distribution around a value that is three lower (log2) than the lowest intensity value recorded for the precursor/peptide and that has a spread of the mean standard deviation for the precursor/peptide.
• "MAR" data is imputed using the mean and variance of the condition with the missing data.
• NA data is not imputed.

Distributions

Summary
Below we can see the distribution of non-imputed data and data imputed with 4 methods:

1. Shifted Distribution based on the data;
2. Imputation using the Random Forest algorithm;
3. Imputation using the k-Nearest Neighbors algorithm;
4. Imputation using the ludovic method.

For each histogram, we can see that the 24th replicate (blue colour) of the KO cell type has a distribution shift to the right. Our goal is to get this distribution closer to the others, and we can do this using different normalization methods.

Mati Code

Random Forest

k-Nearest Neighbors

Ludovic

Normalization

The 24th replicate of the KO cell type has the distribution shifted to the right and should be more similar to the other replicates. Therefore we tried several normalization methods to achieve this goal. We used the following methods:

1. Z-Score
2. Min-Max
3. Median Scaling
4. MAD (Median Absolute Deviation) Scaling
5. Linear Regression
6. VSN (Variance Stabilization Normalization)
7. EigenMS

Z-Score

Z-score normalization standardizes data by subtracting the mean and dividing by the standard deviation. This technique transforms data into a distribution with a mean of 0 and a standard deviation of 1.

Formula: \(\tilde{y}_{ij} = \frac{y_ij - \bar{y_j}}{\theta_j}\), where:

• \(y_ij\) - value of the LFQ;
• \(\bar{y_j}\) - mean of the LFQ values;
• \(\theta_j\) - standard deviation of the LFQ values.

Summary
The Z-SCORE method brought the distribution closer together, but we lost the original scale of the data and the differences are harder to see after normalization.

DISTRIBUTIONS

Min-Max

Min-max normalization scales data to a specified range (usually [0, 1]) by subtracting the minimum value and dividing by the range of values.

Summary
The MIN-MAX method didn’t succeed in bringing the distribution closer together. In conclusion, we probably shouldn’t use this method.

DISTRIBUTIONS

Median Scaling

Median Scaling, also called the Robust Scalar, transforms x to x’ by subtracting each value of features by the median and dividing it by the interquartile range between the 1st quartile (25th quantile) and the 3rd quartile (75th quantile). The formula:

\[x' = \frac{x-median(x)}{(Q3-Q1)}\]

Summary
The MEDIAN scaling method brought the distributions closer together, but again we lost the original scale of the data and the differences are harder to see after normalization.

DISTRIBUTIONS

MAD Scaling

MAD (Median Absoulte Deviation) Scaling transforms \(x\) to \(x'\) in similar way like Median Scaling, but it substractingthe MAD value from the \(x\). The formula of MAD:

\[\tilde{X} = median(X) \\ MAD = median\Big(|X_i - \tilde{X}|\Big)\]

and the formula of scaling is:

\[x' = \frac{x - MAD}{(Q3-Q1)}\]

Summary
The MAD scaling method brought the distributions closer together, but again we lost the original scale of the data and the differences are harder to see after normalization.

DISTRIBUTIONS

Linear Regression

Normalization by linear regression using a residual of the linear model \(y = ax + b\) to normalize the data. In our project we have two variables: independent and dependent. The independent variable is Batch and the dependent variable is LFQ.

Summary
The linear regression method brought the distributions closer together, but again we lost the original scale of the data and the differences are harder to see after normalization.

DISTRIBUTIONS

VSN

VSN - Variance Stabilization Normalization

The data are returned on a \(glog\) scale to base 2. More precisely, the transformed data are subject to the transformation \(glog_2(f(b)*x+a)+c\), where the function \(glog_2(u)=log_2(u+\sqrt{u*u+1}) = asinh(u)/log(2)\) is called the generalised logarithm, the offset \(a\) and the scaling parameter \(b\) are the fitted model parameters (see references), and \(f(x)=exp(x)\) is a parameter transformation that allows ensuring positivity of the factor in front of x while using an unconstrained optimization over \(b\) [source]. The overall offset \(c\) is computed from the \(b\)’s such that for large \(x\) the transformation approximately corresponds to the \(log_2\) function. This is done separately for each stratum, but with the same value across arrays. More precisely, if the element \(b[s,i]\) of the array \(b\) is the scaling parameter for the \(s\)-th stratum and the \(i\)-th array, then \(c[s]\) is computed as \(log2(2*f(mean(b[,i])))\). The offset \(c\) is inconsequential for all differential expression calculations, but many users like to see the data in a range that they are familiar with.

Summary
The VSN method brought the distributions closer together and we didn’t lose the original scale of the data. The method didn’t lose the differences between KO and WT cells at first glance. In addition, the 24th replicate of KO cell doesn’t stand out only in distribution of imputed data with kNN.

DISTRIBUTIONS

EigenMS

EigenMS is a normalization method developed for high-throughput mass spectrometry (MS) proteomics data. It uses eigenvector-based adjustments to remove unwanted variation, such as batch effects, from the data. The method works by identifying systematic biases in the dataset and applying a correction while preserving biologically relevant signals.

In the EigenMS method, eigentrends are calculated by decomposing the systematic variation in the proteomics data matrix using Singular Value Decomposition (SVD). This is typically expressed as:

\[X = U \Sigma V'\] Where:

• \(X\) represents the original data matrix.
• \(U\) and \(v'\) are orthogonal matrices.
• \(\Sigma\) is a diagonal matrix of singular values.

The decomposition isolates systematic trends (eigentrends) as columns in \(U\), which can then be removed to adjust for unwanted variation.

Summary
The EigenMS method was the one which kept the differences between the replicates and cell types (KO and WT). The distribution are close together and we solve the problem of 24th replicate. In my opinion the EigenMS method gives the best results.

DISTRIBUTIONS

Statistics Metrics

CV

COEFFICIENT OF VARIATION (CV) / RELATIVE STANDARD DEVIATION (RSD)

This is a way to measure how spread out values are in a dataset relative to the mean. A lower RSD/CV indicates better normalization. It is calculated as:

\(CV = \frac{\sigma}{\mu}\)

where:

• \(\sigma\): The standard deviation of dataset
• \(\mu\): The mean of dataset

Summary For the biological data the cut-off level of CV is around 20 (+/-5). After close looking at the tables in the **Summary** section, we can safely say that the EigenMS method had the CV around 20 in all imputation methods. In fact, from this point of view we can distinguish four such methods that had good coefficient of variation.

1. EigenMS
2. MAD Scaling
3. Min-Max Normalization
4. Z-Score Normalization

When we take a closer look at the violin plots, we can see the strange pattern in the results for the kNN imputation method. Namely, many values have a CV close to the 0, which can mean that the standard deviation of the values is also close to the 0, which means that the values are very similar, which shouldn’t be the case.

The violin plots for the other methods look better and show the distribution/density of the CV correctly.

Non-Imputed

Mati Code

Random Forest

K-Nearest Neighbors

Ludovic

Summary

Non-Imputed

Imputed - Mati Code

Random Forest

K-Nearest Neighbors

Ludovic

ICC

After normalization, biological replicates should group more tightly. We can assess this by measuring the Intraclass Correlation Coefficient (ICC) to see if replicates cluster together.

A guidelines for interpretation by Koo and Li (2016):

• below 0.50: poor ■
• between 0.50 and 0.75: moderate ■

• between 0.75 and 0.90: good ■
• above 0.90: excellent ■

Summary We can clearly see that the replicates are the most tight after the normalization using the EigenMS method. Almost all methods group replicates more tightly, except the Min-Max and Median methods, where the ICCs were moderate or good.

Non-Imputed

Mati Code

Random Forest

K-Nearest Neighbors

Ludovic

t-test

In our project we used the paired t-test, which is a test based on the differences between the values of a single pair, i.e., one deducted from the other. The formula for the paired t-test is given by:

\[t = \frac{\sum{d}}{\sqrt{\frac{n(\sum{d^2})-(\sum{d})^2}{n-1}}}\]

Where:

• \(d\) is a difference per paired value;
• \(n\) is number of samples.

In our calculations, we tested the following hypothesis:

\(H_0\): There is no significant difference between the means of the two groups.
\(H_1\): There is a significant difference between the two population means.

Raw Data

Summary
The p-value plots show us that we had to normalize the data, because the output below doesn’t show strong differences between KO and WT cells.

Z-Score

Summary
The p-value plots look very bad, which means that we have lost whole differences between the samples.

Min-Max

Summary
Again the p-value plots look very bad, which means that we have again lost whole differences between the samples.

Median Scaling

Summary
Again the p-value plots look very bad, which means that we have again lost whole differences between the samples.

Mad Scaling

Summary
Again the p-value plots look very bad, which means that we have again lost whole differences between the samples.

Linear Regression

Summary
Again the p-value plots look very bad, which means that we have again lost whole differences between the samples.

VSN

Summary
Again the p-value plots look very bad, which means that we have again lost whole differences between the samples.

EigenMS

Summary
Finally, we have the expected p-value plots, which means that the EigenMS method again has the best results.

Fold Change

Raw Data

Non-Imputed

Mati Code

Random Forest

k-Nearest Neighbors

Ludovic

EigenMS

Non-Imputed

Mati Code

Random Forest

k-Nearest Neighbors

Ludovic

Conclusion

Summarizing everything above, we can safely say that the EigenMS method gives the best results. The distributions of the two normalization methods (VSN and EigenMS) were the best until the very end, but in the statistical metrics the EigenMS method won. It is worth noting that the correct shape of the p-value plot was only obtained for the data normalized with the EigenMS method.

About the imputation method.

1. We can safely say that the Random Forest algorithm and the Ludovic method had the best performance.
2. The kNN method corrects the shift in the distribution of the KO 24 replicates, but in the next calculations we discover a strange value in the imputed data. In fact, in some rows the values were the same, making it difficult to calculate the t-test and introducing an error in the CV calculations.
3. The method that samples values from the distribution calculated on the raw data creates the additional peak of columns that doesn’t appear in the raw non-imputed data.

---

References

1. A systematic evaluation of normalization methods in quantitative label-free proteomics
2. How to Normalize Data in R for my Data: Methods and Examples
3. EigenMS
4. Normalization of peak intensities in bottom-up MS-based proteomics using singular value decomposition
5. Metabolomics Data Normalization with EigenMS